%     VIRTUAL SORIC - программа для получения эффективных свойств заполнителей SORIC
%     Copyright (C) 2020  Михаил Козлов
%     Сайт проекта: <http://www.fastenerdoctor.ru/>
%     Адрес электронной почты для связи с разработчиком:
%     <doctor.fastener@gmail.com>
% 
%     Эта программа является свободным программным обеспечением: Вы можете
%     распространять ее и/или изменять, соблюдая условия Универсальной общественной
%     лицензии GNU, опубликованной Фондом свободного программного обеспечения,
%     либо редакции 3, либо (на Ваше усмотрение) любой редакции, выпущенной позже.
% 
%     Эта программа распространяется в надежде на то, что она окажется полезной,
%     но БЕЗ КАКИХ-ЛИБО ГАРАНТИЙ, включая даже подразумеваемую гарантию ТОВАРНОЙ ПРИГОДНОСТИ,
%     либо ПРИГОДНОСТИ ДЛЯ ОПРЕДЕЛЕННЫХ ЦЕЛЕЙ. Ознакомьтесь с Универсальной общественной
%     лицензией GNU для получения более подробной информации.
% 
%     Вы должны были получить копию Универсальной общественной лицензии GNU
%     вместе с этой программой. Если Вы ее не получили, то перейдите по адресу:
%     <http://www.gnu.org/licenses/>.
% 
%     This program is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
% 
%     This program is distributed in the hope that it will be useful,
%     but WITHOUT ANY WARRANTY; without even the implied warranty of
%     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
%     GNU General Public License for more details.
% 
%     You should have received a copy of the GNU General Public License
%     along with this program.  If not, see <https://www.gnu.org/licenses/>.

clear all; %#ok<CLALL>
clc;
more off;

% #######################################################################################################
% ##                                                                                                   ##
% ##                                        НАСТРОЙКИ ПРОГРАММЫ                                        ##
% ##                                                                                                   ##
% #######################################################################################################

% Для включения возможности вывода информации на русском языке при выполнении программы в среде GNU Octave
% необходимо добавить в файл C:\Octave\Octave-4.4.0\share\octave\site\m\startup\octaverc  следующую строку:
% system("chcp 1251",1);

% ----------------------------------- Параметры ячейки периодичности ------------------------------------

% Модель ячейки
Soric_Model                = 'CONSTANT_WEB';    % CONSTANT_WEB | RIGHT_HEX    % Для SORIC LRC следует использовать модель CONSTANT_WEB

% Геометрия ячейки периодичности в плоскости
RVE_Length_1               = 28.19;                % Размер по направлению 1 (направление с параллельными перемычками в листе материала)
RVE_Length_2               = 19.28272727;          % Размер по направлению 2
Web_Volume_Fraction        = 0.025056068;          % Объемная доля перемычек в листе материала SORIC

% Высота ячейки периодичности
RVE_Thickness              = 2.648;                % Толщина листа материала

% Параметры оболочек шестиугольной ячейки (поверхностные слои связующего)
Shell_Layer                = 'TRUE';               % TRUE | FALSE    % Наличие в модели верхней и нижней оболочек шестиугольной ячейки, состоящих из связующего

        % Если Shell_Layer == TRUE
        Shell_Layer_Ratio          = 0.075;                % Относительная толщина одной оболочки
        Wrinkle_Coeff              = 1.5285;               % Коэффициент волнистости оболочек (коэффициент увеличения толщины оболочек при вычислении объемной доли связующего внутри шестиугольной ячейки)

% Объемная доля связующего в шестиугольной ячейке (вместе с оболочками при их наличии)
V_resin_cell               = 0.295012295;

% Суммарная объемная доля связующего во всей ячейке периодичности (для осреднения плотности и теплоемкости)
V_resin_total              = 0.339571163353546;

% ------------------------------------------- Параметры сетки -------------------------------------------

Element_Order              = 'LINEAR'; % LINEAR | QUADRATIC    % Порядок конечных элементов
N_Elements_Edge            = 15; % Число конечных элементов по поперечному ребру шестиугольной ячейки
N_Elements_High            = 8; % Число конечных элементов по высоте ячейки периодичности (без учета оболочек)

% -------------------------------- Параметры массива ячеек периодичности --------------------------------

RVE_Array_1                = 1; % Число ячеек по направлению 1
RVE_Array_2                = 1; % Число ячеек по направлению 2
RVE_Array_3                = 1; % Число ячеек по направлению 3

% ----------------------------------------- Параметры расчета -------------------------------------------

Analysis_Type              = 'MECHANICAL'; % MECHANICAL | HEAT_TRANSFER | MESH_GENERATOR

                           % MECHANICAL - получение эффективных механических свойств с помощью решателя Abaqus
                           % HEAT_TRANSFER - получение эффективных свойств теплопереноса с помощью решателя Abaqus
                           % MESH_GENERATOR - только генерация сетки (без расчета)
                           
Mechanical_BC_TYPE         = 'HARD'; % SOFT | HARD    % Тип граничных условий в механическом расчете

                           % SOFT - периодические граничные условия со свободной границей ячейки
                           % HARD - периодические граничные условия с фиксированной границей ячейки (использовать с большой осторожностью!)
                           
Constraint_Directions      = 'Z';    % Направления, задающие грани с фиксированной границей в случае Mechanical_BC_TYPE == HARD. Строка, состоящая из символов X, Y, Z.

Chemical_Shrinkage         = 'FALSE'; % TRUE | FALSE    % Расчет эффективных коэффициентов химической усадки помимо основных механических свойств

% ----------------------------------------- Свойства материала ------------------------------------------

% ДЛЯ МЕХАНИЧЕСКОГО И ТЕПЛОВОГО РАСЧЕТА

Resin_rho                  = 1.1423e-006;    % Плотность связующего
Filler_rho                 = 0.132664073e-006;    % Плотность сухого материала

% ДЛЯ МЕХАНИЧЕСКОГО РАСЧЕТА

Resin_E                    = 3388.11080952381;    % Модуль упругости связующего
Resin_nu                   = 0.369712485714286;    % Коэффициент Пуассона связующего
Resin_alpha                = 73.4255;    % Коэффициент линейного теплового расширения связующего
Resin_beta                 = 0.01;    % Коэффициент химической усадки связующего

Filler_E                   = 15;    % Модуль упругости сухого материала
Filler_nu                  = 0.37;    % Коэффициент Пуассона сухого материала
Filler_alpha               = 60;    % Коэффициент линейного теплового расширения сухого материала
Filler_beta                = 0;    % Коэффициент химической усадки сухого материала

% ДЛЯ ТЕПЛОВОГО РАСЧЕТА

Resin_Cp                   = 1289.4;    % Удельная теплоемкость связующего
Resin_lambda               = 0.154611954;    % Коэффициент теплопроводности связующего

Filler_Cp                  = 4481;    % Удельная теплоемкость сухого материала
Filler_lambda              = 0.0701341972249683;    % Коэффициент теплопроводности сухого материала

% ----------------------------------------- Параметры запуска -------------------------------------------

Work_Directory             = 'M:\Abaqus_Work';    % Рабочая папка для проведения расчета в Abaqus
CPU                        = 4;                   % Число используемых ядер для расчета в Abaqus

% -------------------------------------------------------------------------------------------------------

% #######################################################################################################
% ##                                                                                                   ##
% ##                                     КОНЕЦ НАСТРОЕК ПРОГРАММЫ                                      ##
% ##                                                                                                   ##
% #######################################################################################################

global ErrCode;
global Merge_Tolerance_Multiplier;
global Element_Critical_Ratio;
global software_version;
global software_date;
global software_year;

ErrCode = 0;
Merge_Tolerance_Multiplier = [];
Element_Critical_Ratio = [];
software_version = [];
software_date = [];
software_year = [];

% Отсчет времени выполнения программы
tic;

% Инициализация базовых переменных
Init_Variables();

% Печать приветствия
Print_Greeting();

% Проверка параметров
switch Soric_Model
    case {'CONSTANT_WEB', 'RIGHT_HEX'}
        % Допустимые значения
    otherwise
        ErrCode = -1;
        fprintf (1, 'Недопустимое значение параметра Soric_Model !\n');
        fprintf (1, 'Разрешенные значения: CONSTANT_WEB, RIGHT_HEX\n');
        fprintf (1, 'Работа программы завершена с ошибкой!\n\n');
        toc;
        return;
end

if (RVE_Length_1 <= 0)
    ErrCode = -1;
    fprintf (1, 'Недопустимое значение параметра RVE_Length_1 !\n');
    fprintf (1, 'Работа программы завершена с ошибкой!\n\n');
    toc;
    return;
end

if (RVE_Length_2 <= 0)
    ErrCode = -1;
    fprintf (1, 'Недопустимое значение параметра RVE_Length_2 !\n');
    fprintf (1, 'Работа программы завершена с ошибкой!\n\n');
    toc;
    return;
end

if ((Web_Volume_Fraction <= 0) || (Web_Volume_Fraction >= 1))
    ErrCode = -1;
    fprintf (1, 'Недопустимое значение параметра Web_Volume_Fraction !\n');
    fprintf (1, 'Работа программы завершена с ошибкой!\n\n');
    toc;
    return;
end

if (RVE_Thickness <= 0)
    ErrCode = -1;
    fprintf (1, 'Недопустимое значение параметра Web_Volume_Fraction !\n');
    fprintf (1, 'Работа программы завершена с ошибкой!\n\n');
    toc;
    return;
end

switch Shell_Layer
    case {'TRUE', 'FALSE'}
        % Допустимые значения
    otherwise
        ErrCode = -1;
        fprintf (1, 'Недопустимое значение параметра Shell_Layer !\n');
        fprintf (1, 'Разрешенные значения: TRUE, FALSE\n');
        fprintf (1, 'Работа программы завершена с ошибкой!\n\n');
        toc;
        return;
end

if ((Shell_Layer_Ratio <= 0) || (Shell_Layer_Ratio >= 0.5))
    ErrCode = -1;
    fprintf (1, 'Недопустимое значение параметра Shell_Layer_Ratio !\n');
    fprintf (1, 'Работа программы завершена с ошибкой!\n\n');
    toc;
    return;
end

if (Wrinkle_Coeff <= 0)
    ErrCode = -1;
    fprintf (1, 'Недопустимое значение параметра Wrinkle_Coeff !\n');
    fprintf (1, 'Работа программы завершена с ошибкой!\n\n');
    toc;
    return;
end

if ((V_resin_cell <= 0) || (V_resin_cell >= 1))
    ErrCode = -1;
    fprintf (1, 'Недопустимое значение параметра V_resin_cell !\n');
    fprintf (1, 'Работа программы завершена с ошибкой!\n\n');
    toc;
    return;
end

if ((V_resin_total <= 0) || (V_resin_total >= 1))
    ErrCode = -1;
    fprintf (1, 'Недопустимое значение параметра V_resin_total !\n');
    fprintf (1, 'Работа программы завершена с ошибкой!\n\n');
    toc;
    return;
end

switch Element_Order
    case {'LINEAR', 'QUADRATIC'}
        % Допустимые значения
    otherwise
        ErrCode = -1;
        fprintf (1, 'Недопустимое значение параметра Element_Order !\n');
        fprintf (1, 'Разрешенные значения: LINEAR, QUADRATIC\n');
        fprintf (1, 'Работа программы завершена с ошибкой!\n\n');
        toc;
        return;
end

if (N_Elements_Edge ~= floor(N_Elements_Edge))
    ErrCode = -1;
    fprintf (1, 'Недопустимое значение параметра N_Elements_Edge !\n');
    fprintf (1, 'Параметр должен быть целым числом.\n');
    fprintf (1, 'Работа программы завершена с ошибкой!\n\n');
    toc;
    return;
end

if (N_Elements_Edge < 1)
    ErrCode = -1;
    fprintf (1, 'Недопустимое значение параметра N_Elements_Edge !\n');
    fprintf (1, 'Параметр должен быть больше или равен единице.\n');
    fprintf (1, 'Работа программы завершена с ошибкой!\n\n');
    toc;
    return;
end

if (N_Elements_High ~= floor(N_Elements_High))
    ErrCode = -1;
    fprintf (1, 'Недопустимое значение параметра N_Elements_High !\n');
    fprintf (1, 'Параметр должен быть целым числом.\n');
    fprintf (1, 'Работа программы завершена с ошибкой!\n\n');
    toc;
    return;
end

if (N_Elements_High < 1)
    ErrCode = -1;
    fprintf (1, 'Недопустимое значение параметра N_Elements_High !\n');
    fprintf (1, 'Параметр должен быть больше или равен единице.\n');
    fprintf (1, 'Работа программы завершена с ошибкой!\n\n');
    toc;
    return;
end

if (RVE_Array_1 ~= floor(RVE_Array_1))
    ErrCode = -1;
    fprintf (1, 'Недопустимое значение параметра RVE_Array_1 !\n');
    fprintf (1, 'Параметр должен быть целым числом.\n');
    fprintf (1, 'Работа программы завершена с ошибкой!\n\n');
    toc;
    return;
end

if (RVE_Array_1 < 1)
    ErrCode = -1;
    fprintf (1, 'Недопустимое значение параметра RVE_Array_1 !\n');
    fprintf (1, 'Параметр должен быть больше или равен единице.\n');
    fprintf (1, 'Работа программы завершена с ошибкой!\n\n');
    toc;
    return;
end

if (RVE_Array_2 ~= floor(RVE_Array_2))
    ErrCode = -1;
    fprintf (1, 'Недопустимое значение параметра RVE_Array_2 !\n');
    fprintf (1, 'Параметр должен быть целым числом.\n');
    fprintf (1, 'Работа программы завершена с ошибкой!\n\n');
    toc;
    return;
end

if (RVE_Array_2 < 1)
    ErrCode = -1;
    fprintf (1, 'Недопустимое значение параметра RVE_Array_2 !\n');
    fprintf (1, 'Параметр должен быть больше или равен единице.\n');
    fprintf (1, 'Работа программы завершена с ошибкой!\n\n');
    toc;
    return;
end

if (RVE_Array_3 ~= floor(RVE_Array_3))
    ErrCode = -1;
    fprintf (1, 'Недопустимое значение параметра RVE_Array_3 !\n');
    fprintf (1, 'Параметр должен быть целым числом.\n');
    fprintf (1, 'Работа программы завершена с ошибкой!\n\n');
    toc;
    return;
end

if (RVE_Array_3 < 1)
    ErrCode = -1;
    fprintf (1, 'Недопустимое значение параметра RVE_Array_3 !\n');
    fprintf (1, 'Параметр должен быть больше или равен единице.\n');
    fprintf (1, 'Работа программы завершена с ошибкой!\n\n');
    toc;
    return;
end

switch Analysis_Type
    case {'MECHANICAL', 'HEAT_TRANSFER', 'MESH_GENERATOR'}
        % Допустимые значения
    otherwise
        ErrCode = -1;
        fprintf (1, 'Недопустимое значение параметра Analysis_Type !\n');
        fprintf (1, 'Разрешенные значения: MECHANICAL, HEAT_TRANSFER, MESH_GENERATOR\n');
        fprintf (1, 'Работа программы завершена с ошибкой!\n\n');
        toc;
        return;
end

switch Mechanical_BC_TYPE
    case {'SOFT', 'HARD'}
        % Допустимые значения
    otherwise
        ErrCode = -1;
        fprintf (1, 'Недопустимое значение параметра Mechanical_BC_TYPE !\n');
        fprintf (1, 'Разрешенные значения: SOFT, HARD\n');
        fprintf (1, 'Работа программы завершена с ошибкой!\n\n');
        toc;
        return;
end

if ((numel(Constraint_Directions) < 1) || (numel(Constraint_Directions) > 3))
    ErrCode = -1;
    fprintf (1, 'Недопустимое значение параметра Constraint_Directions !\n');
    fprintf (1, 'Разрешенные значения: строка, состоящая из символов X, Y, Z.\n');
    fprintf (1, 'Работа программы завершена с ошибкой!\n\n');
    toc;
    return;
end

for i=1:numel(Constraint_Directions)
    switch Constraint_Directions(i)
        case {'X', 'Y', 'Z'}
            % Допустимые значения
        otherwise
            ErrCode = -1;
            fprintf (1, 'Недопустимое значение параметра Constraint_Directions !\n');
            fprintf (1, 'Разрешенные значения: строка, состоящая из символов X, Y, Z.\n');
            fprintf (1, 'Работа программы завершена с ошибкой!\n\n');
            toc;
            return;
    end
end

switch Chemical_Shrinkage
    case {'TRUE', 'FALSE'}
        % Допустимые значения
    otherwise
        ErrCode = -1;
        fprintf (1, 'Недопустимое значение параметра Chemical_Shrinkage !\n');
        fprintf (1, 'Разрешенные значения: TRUE, FALSE\n');
        fprintf (1, 'Работа программы завершена с ошибкой!\n\n');
        toc;
        return;
end

if (CPU ~= floor(CPU))
    ErrCode = -1;
    fprintf (1, 'Недопустимое значение параметра CPU !\n');
    fprintf (1, 'Параметр должен быть целым числом.\n');
    fprintf (1, 'Работа программы завершена с ошибкой!\n\n');
    toc;
    return;
end

if (CPU < 1)
    ErrCode = -1;
    fprintf (1, 'Недопустимое значение параметра CPU !\n');
    fprintf (1, 'Параметр должен быть больше или равен единице.\n');
    fprintf (1, 'Работа программы завершена с ошибкой!\n\n');
    toc;
    return;
end

% Счетчики числа сообщений
N_Info = 0;
N_Warnings = 0;

% Анализ настроек программы
if (strcmp(Analysis_Type, 'MECHANICAL'))
    
    if (strcmp(Mechanical_BC_TYPE, 'HARD'))
        fprintf (1, 'ПРЕДУПРЕЖДЕНИЕ! В расчете используются периодические граничные условия с фиксированной границей\n');
        fprintf (1, '                ячейки (Mechanical_BC_TYPE = HARD). Требуется особое внимание при анализе\n');
        fprintf (1, '                полученных эффективных свойств.\n\n');
        N_Warnings = N_Warnings + 1;
    end
    
end

% Толщина поверхностного слоя связующего
Shell_Layer_Thickness = [];

if (strcmp(Shell_Layer, 'TRUE'))
    Shell_Layer_Thickness = Shell_Layer_Ratio*RVE_Thickness;
end

% Объемная доля связующего внутри шестиугольной ячейки
if (strcmp(Shell_Layer, 'TRUE'))
    V_resin_cell = (V_resin_cell - 2*Shell_Layer_Ratio*Wrinkle_Coeff)/(1 - 2*Shell_Layer_Ratio*Wrinkle_Coeff);
end

% Создание базовых точек
Points = Construct_Geometry(Soric_Model, RVE_Length_1, RVE_Length_2, Web_Volume_Fraction);

if (ErrCode < 0)
    fprintf (1, 'Работа программы завершена с ошибкой!\n\n');
    toc;
    return;
end

% Размеры одинарной ячейки периодичности
RVE_Size = zeros(1,3);
RVE_Size(1) = Points(23,1) - Points(29,1);
RVE_Size(2) = Points(31,2) - Points(26,2);
RVE_Size(3) = RVE_Thickness;

% Построение сетки
[NODES, ELEMENTS, RVE_Coordinates, Node_Merge_Tolerance] = Generate_Mesh(Points, N_Elements_Edge, RVE_Thickness, N_Elements_High, Shell_Layer, Shell_Layer_Thickness, RVE_Size, Element_Order, RVE_Array_1, RVE_Array_2, RVE_Array_3);

if (ErrCode < 0)
    fprintf (1, 'Работа программы завершена с ошибкой!\n\n');
    toc;
    return;
end

% Корректировка Work_Directory
if (Work_Directory(length(Work_Directory)) == '\')
    Work_Directory = Work_Directory(1:length(Work_Directory)-1);
end

% Заполнение структуры со свойствами материала
Material_Properties.Resin_E = Resin_E;
Material_Properties.Resin_nu = Resin_nu;
Material_Properties.Resin_alpha = Resin_alpha;
Material_Properties.Resin_beta = Resin_beta;
Material_Properties.Resin_rho = Resin_rho;
Material_Properties.Resin_Cp = Resin_Cp;
Material_Properties.Resin_lambda = Resin_lambda;

Material_Properties.Filler_E = Filler_E;
Material_Properties.Filler_nu = Filler_nu;
Material_Properties.Filler_alpha = Filler_alpha;
Material_Properties.Filler_beta = Filler_beta;
Material_Properties.Filler_rho = Filler_rho;
Material_Properties.Filler_Cp = Filler_Cp;
Material_Properties.Filler_lambda = Filler_lambda;

Material_Properties.Cell_E = [];
Material_Properties.Cell_nu = [];
Material_Properties.Cell_alpha = [];
Material_Properties.Cell_beta = [];
Material_Properties.Cell_rho = [];
Material_Properties.Cell_Cp = [];
Material_Properties.Cell_lambda = [];

% Проверка свойств материала
Check_Material_Properties(Material_Properties);

if (ErrCode < 0)
    fprintf (1, 'Работа программы завершена с ошибкой!\n\n');
    toc;
    return;
end

% Ветвление по типу анализа
switch Analysis_Type
    
    case 'MECHANICAL'
        
        Material_Properties = Mechanical_Procedure(Work_Directory, CPU, NODES, ELEMENTS, Material_Properties, RVE_Coordinates, Node_Merge_Tolerance, Chemical_Shrinkage, Mechanical_BC_TYPE, V_resin_cell, V_resin_total, Shell_Layer, Constraint_Directions);
        
        if (ErrCode < 0)
            fprintf (1, 'Работа программы завершена с ошибкой!\n\n');
            toc;
            return;
        end
        
    case 'HEAT_TRANSFER'
        
        Material_Properties = Heat_Transfer_Procedure(Work_Directory, CPU, NODES, ELEMENTS, Material_Properties, RVE_Coordinates, Node_Merge_Tolerance, V_resin_cell, V_resin_total, Shell_Layer);
        
        if (ErrCode < 0)
            fprintf (1, 'Работа программы завершена с ошибкой!\n\n');
            toc;
            return;
        end
        
    case 'MESH_GENERATOR'
        
        Material_Properties = Compute_Cell_Mech_Props(Material_Properties, V_resin_cell);
        
        if (ErrCode < 0)
            fprintf (1, 'Работа программы завершена с ошибкой!\n\n');
            toc;
            return;
        end
        
        Write_Input_Mesh(Work_Directory, NODES, ELEMENTS, Material_Properties, Shell_Layer);

        if (ErrCode < 0)
            fprintf (1, 'Работа программы завершена с ошибкой!\n\n');
            toc;
            return;
        end
    
end

% Печать разделителя
Print_Delimiter();

% Сообщение об успешном завершении работы программы
fprintf (1, 'Работа программы завершена успешно!\n\n');

% Вывод числа предупреждений
if (N_Warnings > 0)
    fprintf (1, '\tПолучено предупреждений: %i\n', N_Warnings);
end

% Вывод числа информационных сообщений
if (N_Info > 0)
    fprintf (1, '\tПолучено информационных сообщений: %i\n', N_Info);
end

% Печать новой строки
if ((N_Warnings > 0) || (N_Info > 0))
    fprintf (1, '\n');
end

% Вывод времени выполнения программы
toc;
